! This is a 1D advection example using square initial condition and periodic
! boundary condition for upwind finite difference scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-20: Initial creation.
! - 2019-03-13: Add parameter beta to control the upwind degree.
! - 2024-03-21: Use adv_1d_test_case_mod.

program upwind_adv_1d_case

  use adv_1d_test_case_mod

  implicit none

  real, allocatable :: rho(:,:)     ! Tracer density being advected at cell centers
  real, allocatable :: flx(:)       ! Flux at cell interfaces
  real, parameter :: beta = 1.0     ! Upwind weight
  integer, parameter :: ns = 1      ! Stencil width
  integer i
  character(30), parameter :: scheme = 'upwind'

  namelist /params/ nx, nt, dt, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,2))
  allocate(flx(1-ns:nx+ns))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, 0, ns, nx, x, rho(:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call upwind(rho(:,old), flx)
    do i = 1, nx
      rho(i,new) = rho(i,old) - dt / dx * (flx(i) - flx(i-1))
    end do
    call apply_bc(ns, nx, rho(:,new))
    ! Change time indices.
    i = old; old = new; new = i
    time_step = time_step + 1
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    print *, time_step, sum(rho(1:nx,old))
  end do

  call adv_1d_test_case_final()

  deallocate(rho)
  deallocate(flx)

contains

  subroutine upwind(q, f)

    real, intent(in ) :: q(1-ns:nx+ns)
    real, intent(out) :: f(1-ns:nx+ns)

    integer i

    do i = 1, nx
      f(i) = 0.5d0 * (u * (q(i+1) + q(i)) - beta * abs(u) * (q(i+1) - q(i)))
    end do
    call apply_bc(ns, nx, flx)

  end subroutine upwind

end program upwind_adv_1d_case
